library(ggplot2)
library(tidyverse)
library(rdd)
library(rddensity)
library(rdrobust)
library(rddtools)
library(xtable)
library(stargazer)



########################################################################################################
## Read in ballot measures data

taxes <- read.csv("data/taxes.csv")


# 448 cities voted on at least one tax increase over the panel
length(unique(taxes$city))


# How many elections per city that where tax measures were on the ballot?
# Average: 4 (Max: 34 in San Francisco)
tax.elections <- taxes %>% group_by(city) %>% summarize(n = n())
summary(tax.elections$n)
tax.elections[tax.elections$n==34,]



########################################################################################################
## Figure 1: Histogram of Proposed Local Tax Measures

taxes <- read.csv("data/taxes.csv")

## Re-classify uncommon taxes as Misc for plotting
plot.dta <- taxes
plot.dta$tax <- ifelse(plot.dta$reclassification=="Gasoline Tax" | 
                         plot.dta$reclassification=="Development Tax", 
                       "Miscellaneous Tax", plot.dta$reclassification)


(p <- ggplot(plot.dta, aes(x = factor(tax))) +
    geom_histogram(stat = "count", width = 0.5, fill = "steelblue", color = "black") +
    labs(y = "Count", x = "") +
    theme_minimal() +  
    theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
          text=element_text(size=12, family="serif")))

ggsave("output/figure1.pdf", plot = p,
       width = 5, height = 3.5, units = "in", device = "pdf")


########################################################################################################
## Figure 2: Tax Measure Passage by Year

taxes <- read.csv("data/taxes.csv")

plot.dta <- taxes
plot.dta$tax <- ifelse(plot.dta$reclassification=="Gasoline Tax" | 
                         plot.dta$reclassification=="Development Tax", 
                       "Miscellaneous Tax", plot.dta$reclassification)

plot.dta$outcome <- factor(plot.dta$outcome, levels = c(2,1), labels = c("Fail", "Pass"))
plot.dta <- plot.dta[!is.na(plot.dta$outcome),]

(p <- ggplot(plot.dta, aes(x = factor(year), fill = outcome)) +
    geom_bar(stat="count", color = "black") + scale_x_discrete(breaks = seq(1996, 2022, 2)) +
    labs(x = "Year", y = "Count", fill = "Outcome") +
    scale_fill_brewer(palette = "Paired") +
    theme_bw(base_family = "serif", base_size = 12) +
    theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
          legend.position="top"))


ggsave("output/figure2.pdf", plot = p,
       width = 5.5, height = 3.5, units = "in", device = "pdf")


########################################################################################################
## Figure 3a: McCrary (2008) sorting test 

taxes <- read.csv("data/taxes.csv")

pdf(file = "output/figure3a.pdf",
    width = 4.7, height = 3.8)

par(mar = c(5, 5, 2, 2), family = "serif", tck = -.018, las = 1, mgp = c(3, 0.6, 0))
DCdensity(taxes$marg, c = 0)
title(cex.lab = 1.2, xlab = "Percent Supporting Tax Increase", ylab = "Density") 
abline(v = 0, lty = 2)

dev.off()

## Figure 3b: Cattaneo, Jansson, and Ma (2020) sorting test 

denstest_1 <- rddensity(taxes$marg, c = 0, p=3)
summary(denstest_1)

plot <- rdplotdensity(denstest_1, taxes$marg, xlabel = "\n Percent Supporting Tax Increase",
                      lcol=1, histFillCol=1, histBreaks = seq(-.4, .4, .02), ylabel = "Density") 
(p <- plot$Estplot + geom_vline(xintercept = 0, linetype = 2) + xlim(-.4, .4) +
    theme_bw(base_family = "serif", base_size = 14) + theme(legend.position = "none"))

ggsave("output/figure3b.pdf", plot = p,
       width = 4.5, height = 3.8, units = "in", device = "pdf")


########################################################################################################
## Figure 4a: Tax Increases Do Not Change Electoral Competition

dta <- read.csv("data/rdd-data.csv")

plot <- rdplot(dta$prop.running.change, dta$marg, binselect = "qs", # x.lim = c(-.3, .3),
               #  y.lim = c(13, 20), x.lim = c(-.4, .4),
               x.label = "Vote Share Margin for Tax Increase", 
               y.label = "Change in Prop. Incumbents Running for Office", 
               col.lines = "black", title = "")

(p <- plot$rdplot + geom_point(aes(x=dta$marg, y = dta$prop.running.change), alpha = .3) +
    theme_bw(base_family = "serif", base_size = 12))

ggsave("output/figure4a.pdf", plot = p,
       width = 4, height = 4, units = "in", device = "pdf")


########################################################################################################
## Figure 4b: Tax Increases Do Not Change Electoral Competition

dta <- read.csv("data/rdd-data.csv")

plot <- rdplot(dta$challengers.change, dta$marg, binselect = "qs", # x.lim = c(-.3, .3),
               #    y.lim = c(-10, 10), #x.lim = c(-.4, .4),
               x.label = "Vote Share Margin for Tax Increase", 
               y.label = "Change in Number of Challengers Per Seat", 
               col.lines = "black", title = "")

(p <- plot$rdplot + geom_point(aes(x=dta$marg, y = dta$challengers.change), alpha = .3) +
    theme_bw(base_family = "serif", base_size = 12))

ggsave("output/figure4b.pdf", plot = p,
       width = 4, height = 4, units = "in", device = "pdf")


########################################################################################################
## Table 1: Minimal Electoral Punishment for Incumbents After Tax Increase

dta <- read.csv("data/rdd-data.csv")

m <- (rdrobust(dta$avg.inc.voteshare, dta$marg, c=0, all=T, p=1,
               cluster = dta$city))

m1 <- (rdrobust(dta$avg.inc.voteshare, dta$marg, c=0, all=T, p=2,
                cluster = dta$city))

m2 <- (rdrobust(dta$avg.inc.voteshare, dta$marg, c=0, all=T, p=3,
                cluster = dta$city))

m3 <- (rdrobust(dta$prop.inc.win, dta$marg, c=0, all=T, p=1,
                cluster = dta$city))

m4 <- (rdrobust(dta$prop.inc.win, dta$marg, c=0, all=T, p=2,
                cluster = dta$city))

m5 <- (rdrobust(dta$prop.inc.win, dta$marg, c=0, all=T, p=3,
                cluster = dta$city))

coef <- round(c(m$coef[2], m1$coef[2], m2$coef[2], m3$coef[2], m4$coef[2], m5$coef[2]), 3)

se <- c(paste0("(",round(m$se[2], 3),")"), 
        paste0("(",round(m1$se[2], 3),")"), 
        paste0("(",round(m2$se[2], 3),")"), 
        paste0("(",round(m3$se[2], 3),")"), 
        paste0("(",round(m4$se[2], 3),")"), 
        paste0("(",round(m5$se[2], 3),")"))

p <- round(c(m$pv[2], m1$pv[2], m2$pv[2], m3$pv[2], m4$pv[2], m5$pv[2]), 2)

bandwidth <- round(c(m$bws[1], m1$bws[1], m2$bws[1], m3$bws[1], m4$bws[1], m5$bws[1]), 3)


obs <- c(sum(m$N), sum(m1$N), sum(m2$N), sum(m3$N))

eff.obs <- c(sum(m$N_h), sum(m1$N_h), sum(m2$N_h), sum(m3$N_h))

means <- round(c(mean(dta$avg.inc.voteshare, na.rm = T), 
                 mean(dta$avg.inc.voteshare, na.rm = T), 
                 mean(dta$avg.inc.voteshare, na.rm = T),
                 mean(dta$prop.inc.win, na.rm = T), 
                 mean(dta$prop.inc.win, na.rm = T), 
                 mean(dta$prop.inc.win, na.rm = T)), 2)

poly <- c(1, 2, 3, 1, 2, 3)

tab <- rbind(coef, SE = se, poly = poly, bandwidth = bandwidth, 
             obs = obs, eff.obs = eff.obs, means = means)


columns <- c("(1)", "(2)", "(3)", "(4)", "(5)", "(6)")

colnames(tab) <- columns

rownames(tab) <- c("Tax Increase Approved", "", "Polynomial", "Bandwidth", 
                   "Observations", "Obs. in Bandwidth", "Mean Outcome")

tab <- xtable(tab, align = c("l", "c", "c", "c", "c", "c", "c"))


addtorow <- list()
addtorow$pos <- list(-1,0, 2,nrow(tab))
addtorow$command <- c("\\\\[-1.8ex] \\hline \\hline \\\\[-1.8ex]  & 
                      \\multicolumn{3}{c}{Incumbent Vote Share} & 
                      \\multicolumn{3}{c}{Prop. Incumbents Winning} 
                      \\\\ \\cline{2-4} \\cline{5-7} \\\\[-5ex] \\\\",
                      "\\hline \\\\[-1.8ex] ", "\\\\[-1.8ex] \\hline \\\\[-1.8ex]", 
                      "\\\\[-1.8ex] \\hline \\hline \\\\[-1.8ex]")

print(tab, add.to.row = addtorow, hline.after = NULL,
      floating = FALSE, file = "output/table1.tex")



####################################################################################################
## Figure 5: Effect of Tax Increases on Incumbent Vote Share

dta <- read.csv("data/rdd-data.csv")

m <- (rdrobust(dta$avg.inc.voteshare, dta$marg, c=0, all=T, cluster = dta$city, p = 1))

business <- dta[dta$reclassification =="Business Tax",]
m1 <- (rdrobust(business$avg.inc.voteshare, business$marg, c=0, 
                all=T, cluster = business$city, p = 1))

property <- dta[dta$reclassification=="Property Tax",]

m2 <- (rdrobust(property$avg.inc.voteshare, property$marg, c=0, 
                all=T, cluster = property$city, p = 1))

sales <- dta[dta$reclassification=="Sales Tax",]
m3 <- (rdrobust(sales$avg.inc.voteshare, sales$marg, c=0, all=T, 
                cluster = sales$city, p = 1))

utility <- dta[dta$reclassification=="Utility Tax",]
m4 <- (rdrobust(utility$avg.inc.voteshare, utility$marg, c=0, 
                all=T, cluster = utility$city, p = 1))

hotel <- dta[dta$reclassification=="Transient Occupancy Tax",]
m5 <- (rdrobust(hotel$avg.inc.voteshare, hotel$marg, c=0, all=T, cluster = hotel$city, p = 1))


N <- c(sum(m$N), sum(m1$N), sum(m2$N), sum(m3$N), sum(m4$N), sum(m5$N))

modelcoef <- (c(m$coef[2], m1$coef[2], m2$coef[2], m3$coef[2], m4$coef[2], m5$coef[2]))

modelse <- c(m$se[2], m1$se[2], m2$se[2], m3$se[2], m4$se[2], m5$se[2])


ylo = modelcoef - qt(.975, N)*(modelse)
yhi = modelcoef + qt(.975, N)*(modelse)

names = c("Any Tax Increase", "Business Tax Increase", "Property Tax Increase", 
          "Sales Tax Increase", "Utility Tax Increase", "Transient Occupancy Tax Increase")

dfplot = data.frame(names, modelcoef, modelse, ylo, yhi)


(p <- ggplot(dfplot, aes(names, modelcoef, ymin=ylo, ymax=yhi)) +
    geom_pointrange(size = .25) + 
    coord_flip() + geom_hline(yintercept = 0, lty=2) + xlab("") + ylab("") +
    theme_bw(base_family = "serif", base_size = 13)  +
    theme(plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank()))

ggsave("output/figure5.pdf", plot = p,
       width = 5, height = 3.5, units = "in", device = "pdf")


########################################################################################################
# Table 2: Incumbents Are Punished for Business Tax Increases

dta <- read.csv("data/rdd-data.csv")
business <- dta[dta$reclassification == "Business Tax", ]

m <- (rdrobust(business$avg.inc.voteshare, business$marg, c=0, all=T, p=1,
               cluster = business$city))

m1 <- (rdrobust(business$avg.inc.voteshare, business$marg, c=0, all=T, p=2,
                cluster = business$city))

m2 <- (rdrobust(business$avg.inc.voteshare, business$marg, c=0, all=T, p=3,
                cluster = business$city))

m3 <- (rdrobust(business$prop.inc.win, business$marg, c=0, all=T, p=1,
                cluster = business$city))

m4 <- (rdrobust(business$prop.inc.win, business$marg, c=0, all=T, p=2,
                cluster = business$city))

m5 <- (rdrobust(business$prop.inc.win, business$marg, c=0, all=T, p=3,
                cluster = business$city))

coef <- round(c(m$coef[2], m1$coef[2], m2$coef[2], m3$coef[2], m4$coef[2], m5$coef[2]), 3)

se <- c(paste0("(",round(m$se[2], 3),")"), 
        paste0("(",round(m1$se[2], 3),")"), 
        paste0("(",round(m2$se[2],3),")"), 
        paste0("(",round(m3$se[2],3),")"), 
        paste0("(",round(m4$se[2],3),")"), 
        paste0("(",round(m5$se[2],3),")"))

p <- round(c(m$pv[2], m1$pv[2], m2$pv[2], m3$pv[2], m4$pv[2], m5$pv[2]), 2)

bandwidth <- round(c(m$bws[1], m1$bws[1], m2$bws[1], m3$bws[1], m4$bws[1], m5$bws[1]), 3)


obs <- c(sum(m$N), sum(m1$N), sum(m2$N), sum(m3$N))

eff.obs <- c(sum(m$N_h), sum(m1$N_h), sum(m2$N_h), sum(m3$N_h))

means <- round(c(mean(business$avg.inc.voteshare, na.rm = T), 
                 mean(business$avg.inc.voteshare, na.rm = T), 
                 mean(business$avg.inc.voteshare, na.rm = T),
                 mean(business$prop.inc.win, na.rm = T), 
                 mean(business$prop.inc.win, na.rm = T),
                 mean(business$prop.inc.win, na.rm = T)), 2)

poly <- c(1, 2, 3, 1, 2, 3)


tab <- rbind(coef, SE = se, poly = poly, bandwidth = bandwidth, 
             obs = obs, eff.obs = eff.obs, means = means)


columns <- c("(1)", "(2)", "(3)", "(4)", "(5)", "(6)")

colnames(tab) <- columns

rownames(tab) <- c("Business Tax Approved", "", "Polynomial", "Bandwidth", 
                   "Observations", "Obs. in Bandwidth", "Mean Outcome")

tab <- xtable(tab, align = c("l", "c", "c", "c", "c", "c", "c"))


addtorow <- list()
addtorow$pos <- list(-1,0, 2,nrow(tab))
addtorow$command <- c("\\\\[-1.8ex] \\hline \\hline \\\\[-1.8ex]  & 
                      \\multicolumn{3}{c}{Incumbent Vote Share} & 
                      \\multicolumn{3}{c}{Prop. Incumbents Winning} 
                      \\\\ \\cline{2-4} \\cline{5-7} \\\\[-5ex] \\\\",
                      "\\hline \\\\[-1.8ex] ", "\\\\[-1.8ex] \\hline \\\\[-1.8ex]", 
                      "\\\\[-1.8ex] \\hline \\hline \\\\[-1.8ex]")

print(tab, add.to.row = addtorow, hline.after = NULL,
      floating = FALSE, file = "output/table2.tex")



########################################################################################################
## Figure 6a: Effect of Business Tax Increase Across Bandwidths

dta <- read.csv("data/rdd-data.csv")
business <- dta[dta$reclassification == "Business Tax", ]

random.bandwidth <- seq(0.05, .2, .01)

coef <- NA
se <- NA
p <- NA

i <- 1

for(i in 1:length(random.bandwidth)) {
  m <- (rdrobust(business$change.voteshare, business$marg, c=0, h = random.bandwidth[i],
                 all=T, cluster = business$city, p = 1))
  coef <- append(coef, m$coef[2], after = length(coef))
  se <- append(se, m$se[2], after = length(se))
  p <- append(p, m$pv[2], after = length(p))
}

plot.dta <- na.omit(data.frame(coef = coef, se = se, p = p))
plot.dta$bandwidth <- random.bandwidth


(p <- ggplot(plot.dta, aes(x = bandwidth, y = coef)) +
    geom_line() + xlab("Bandwidth") + ylab("Effect on Average Incumbent Voteshare") +
    geom_ribbon(aes(x = bandwidth, ymin = coef - (se*1.96), ymax = coef + se*1.96), alpha = .3) +
    geom_hline(yintercept = 0, linetype = 2) + geom_vline(xintercept = .005, linetype = 1) +
    xlim(0.05, 0.2) +
    theme_bw(base_size = 12, base_family = "serif") +
    theme(plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank()))

ggsave("output/figure6a.pdf", plot = p,
       width = 4, height = 3.1, units = "in", device = "pdf")



########################################################################################################
## Figure 6b: Effect of Business Tax Increase at Random Cut-offs

dta <- read.csv("data/rdd-data.csv")
business <- dta[dta$reclassification == "Business Tax", ]

random.cut <- seq(-.06, .1, .01)

coef <- NA
se <- NA
p <- NA

i <- 1


for(i in 1:length(random.cut)) {
  m <- (rdrobust(business$change.voteshare, business$marg, random.cut[i], 
                 all=T, cluster = business$city, p = 1))
  
  coef <- append(coef, m$coef[2], after = length(coef))
  se <- append(se, m$se[2], after = length(se))
  p <- append(p, m$pv[2], after = length(p))
}

plot.dta <- na.omit(data.frame(coef = coef, se = se, p = p))
plot.dta$cut <- random.cut

(p <- ggplot(plot.dta, aes(x = cut, y = coef)) +
    geom_point() + xlab("Vote Margin Cutoff") + ylab("Effect on Average Incumbent Voteshare") +
    geom_errorbar(aes(x = cut, ymin = coef - (se*1.96), ymax = coef + se*1.96), alpha = .3) +
    geom_hline(yintercept = 0, linetype = 2) + geom_vline(xintercept = .005, linetype = 1) +
    geom_vline(xintercept = -.005, linetype = 1) +
    theme_bw(base_size = 12, base_family = "serif") +
    theme(plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank()))

ggsave("output/figure6b.pdf", plot = p,
       width = 4, height = 3.1, units = "in", device = "pdf")


########################################################################################################
# Table 3: Effect of Business Tax Punishment by Chamber of Commerce Presence 

dta <- read.csv("data/rdd-data.csv")
chambers <- read.csv("data/chambers-of-commerce.csv")

# Look at distribution of chamber type by city
prop.table(table(chambers$chamber_type))

# Restrict to business taxes
business <- dta[dta$reclassification =="Business Tax",]
business <- merge(business, chambers, by = "city")


## Split sample
strong.business <- business[business$chamber_type == "city", ]
weak.business <- business[business$chamber_type != "city", ]


## Effect exclusively driven by cities with their own local chamber of commerce
m <- (rdrobust(strong.business$avg.inc.voteshare, strong.business$marg, c=0, all=T, p=1,
               cluster = strong.business$city))

m1 <- (rdrobust(strong.business$avg.inc.voteshare, strong.business$marg, c=0, all=T, p=2,
                cluster = strong.business$city))


## No effect in cities without a chamber of commerce
m2 <- (rdrobust(weak.business$avg.inc.voteshare, weak.business$marg, c=0, all=T, p=1,
                cluster = weak.business$city))

m3 <- (rdrobust(weak.business$avg.inc.voteshare, weak.business$marg, c=0, all=T, p=2,
                h = .2, cluster = weak.business$city))



coef <- round(c(m$coef[2], m1$coef[2], m2$coef[2], m3$coef[2]), 3)

se <- c(paste0("(",round(m$se[2], 3),")"), 
        paste0("(",round(m1$se[2], 3),")"),
        paste0("(",round(m2$se[2], 3),")"),
        paste0("(",round(m3$se[2], 3),")"))

p <- round(c(m$pv[2], m1$pv[2], m2$pv[2], m3$pv[2]), 2)

bandwidth <- round(c(m$bws[1], m1$bws[1], m2$bws[1], m3$bws[1]), 3)


obs <- c(sum(m$N), sum(m1$N), sum(m2$N), sum(m3$N))

eff.obs <- c(sum(m$N_h), sum(m1$N_h), sum(m2$N_h), sum(m3$N_h))

means <- round(c(mean(strong.business$avg.inc.voteshare, na.rm = T),
                 mean(strong.business$avg.inc.voteshare, na.rm = T),
                 mean(weak.business$avg.inc.voteshare, na.rm = T),
                 mean(weak.business$avg.inc.voteshare, na.rm = T)), 2)

poly <- c(1, 2, 1, 2)


tab <- rbind(coef, SE = se, poly = poly, bandwidth = bandwidth, obs = obs, 
             eff.obs = eff.obs, means = means)


columns <- c("(1)", "(2)", "(3)", "(4)")

colnames(tab) <- columns

rownames(tab) <- c("Business Tax Approved", "", "Polynomial", "Bandwidth", 
                   "Observations", "Obs. in Bandwidth", "Mean Outcome")

tab <- xtable(tab, align = c("l", "c", "c", "c", "c"))


addtorow <- list()
addtorow$pos <- list(-1,0, 2,nrow(tab))
addtorow$command <- c("\\\\[-1.8ex] \\hline \\hline \\\\[-1.8ex]  & 
                      \\multicolumn{4}{c}{Incumbent Vote Share}
                      \\\\ \\cline{2-3} \\cline{4-5} \\\\[-5ex] \\\\ &
                      \\multicolumn{2}{c}{City Chamber of Commerce} & 
                      \\multicolumn{2}{c}{No Chamber of Commerce} \\\\",
                      "\\hline \\\\[-1.8ex] ", "\\\\[-1.8ex] \\hline \\\\[-1.8ex]", 
                      "\\\\[-1.8ex] \\hline \\hline \\\\[-1.8ex]")

print(tab, add.to.row = addtorow, hline.after = NULL,
      floating = FALSE, file = "output/table3.tex")



########################################################################################################
## Figure 7: Effect of Tax Increases on Business Leaders Running for Office

dta <- read.csv("data/rdd-data.csv")

business <- dta[dta$reclassification=="Business Tax",]
sales <- dta[dta$reclassification=="Sales Tax",]
property <- dta[dta$reclassification=="Property Tax",]
utility <- dta[dta$reclassification=="Utility Tax",]
hotel <- dta[dta$reclassification=="Transient Occupancy Tax",]

m <- (rdrobust(dta$business.owners, dta$marg, c=0, all=T, p=1,
               cluster = dta$city))

m1 <- (rdrobust(business$business.owners, business$marg, c=0, all=T, p=1, 
                cluster = business$city))

m2 <- (rdrobust(property$business.owners, property$marg, c=0, all=T, p=1,
                cluster = property$city))

m3 <- (rdrobust(sales$business.owners, sales$marg, c=0, all=T, p=1,
                cluster = sales$city))

m4 <- (rdrobust(hotel$business.owners, hotel$marg, c=0, all=T, p=1,
                cluster = hotel$city))

m5 <- (rdrobust(utility$business.owners, utility$marg, c=0, all=T, p=1,
                cluster = utility$city))


N <- c(sum(m$N), sum(m1$N), sum(m2$N), sum(m3$N), sum(m4$N), sum(m5$N))

modelcoef <- (c(m$coef[2], m1$coef[2], m2$coef[2], m3$coef[2], m4$coef[2], m5$coef[2]))

modelse <- c(m$se[2], m1$se[2], m2$se[2], m3$se[2], m4$se[2], m5$se[2])


ylo = modelcoef - qt(.975, N)*(modelse)
yhi = modelcoef + qt(.975, N)*(modelse)


names = c("Any Tax Increase", "Business Tax Increase", "Property Tax Increase", 
          "Sales Tax Increase", "Transient Occupancy Tax Increase", "Utility Tax Increase")


dfplot = data.frame(names, modelcoef, modelse, ylo, yhi)

(p <- ggplot(dfplot, aes(names, modelcoef, ymin=ylo, ymax=yhi)) +
    geom_pointrange(size = .25) + 
    coord_flip() + geom_hline(yintercept = 0, lty=2) + xlab("") + ylab("") +
    theme_bw(base_family = "serif", base_size = 13)  +
    theme(plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank()))

ggsave("output/figure7.pdf", plot = p,
       width = 5, height = 3.5, units = "in", device = "pdf")

